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ABSTRACT 

Hydrodynamical simulations of star formation have stimulated a need to develop fast and 
robust algorithms for evaluating radiative cooling. Here we undertake a critical evaluation 
of what is currently a popular method for prescribing cooling in Smoothed Particle Hydro- 
dynamical (SPH) simulations, i.e. the poly tropic cooling due originally to Stamatellos et al. 
This method uses the local density and potential to estimate the column density and optical 
depth to each particle and then uses these quantities to evaluate an approximate expression for 
the net radiative cooling. We evaluate the algorithm by considering both spherical and disc- 
like systems with analytic density and temperature structures. In spherical systems, the total 
cooling rate computed by the method is within around 20% for the astrophy sic ally relevant 
case of opacity dominated by ice grains and is correct to within a factor of order unity for a 
range of opacity laws. In disc geometry, however, the method systematically under-estimates 
the cooling by a large factor at all heights in the disc. For the self-gravitating disc studied, 
we find that the method under-estimates the total cooling rate by a factor of ~ 200. This 
discrepancy may be readily traced to the method's systematic over-estimate of the disc col- 
umn density and optical depth, since (being based only on the local density and potential) it 
does not take into account the low column density route for photon escape normal to the disc 
plane. We note that the discrepancy quoted above applies in the case that the star's potential 
is not included in the column density estimate and that even worse agreement is obtained if 
the full (star plus disc) potential is employed. These results raise an obvious caution about 
the method's use in disc geometry whenever an accurate cooling rate is required, although we 
note that there are situations where the discrepancies highlighted above may not significantly 
affect the global outcome of simulations. Finally, we draw attention to our introduction of an 
analytic self-gravitating disc structure that may be of use in the calibration of future cooling 
algorithms. 

Key words: accretion, accretion discs - circumstellar matter - radiative transfer - planetary 
systems: protoplanetary discs - stars: formation - stars: pre-main sequence 



1 INTRODUCTION 

Star formation simulations require computationally efficient yet ro- 
bust algorithms for predicting the thermal evolution of the gas dur- 
ing collapse and fragmentation. A cheap and widely used approach 
(particularly in SPH simulations) is to simply adopt a barotropic 
equation of state, whose form is motivated by the results of spheri- 
cal collapse calculations that include a treatmen t of radiative trans- 
fer dMasunaga & Inutsuka [2000l ; iLarsonl 12005 ) . This prescription 
correctly captures three broad phases of the thermal evolution of 
the gas, i.e. an isothermal stage at low density followed by succes- 
sive stages where the gas heats up in response to increasing density. 
However, by associating gas of a particular density with a corre- 
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sponding temperature it naturally ignores a variety of effects that 
in practice will affect the thermal state of the gas {e.g. the geome- 
try, the temperature of the surrounding gas and the local radiation 
field). 

At the other end of the spectrum of complexity, there is 
a long history of grid based hydrodynamic simulations that fol- 
low the thermal evolution of the gas by including energy trans- 
fer between the gas and radiation field. Such calculations vary 
in the level of sophistication with which they model the radia- 
tion fiel d, ranging from full frequency dependent radiative trans - 
fer (e.g. IWolfire & Cassinellil dl987l) ; iVorke & Sonnhalterl d2002h ) 
to the widely employed and cheaper exp e dient of grey flux l i mited 
diffusi o n (e.g.lBodenheimer et al.lfl99oh:lBolev etaild2007t) ; iBossI 
d2008l) : ICai et alj d2008l) : iKrumholz etafl d2009l) ). More recently, 
several groups have implemented the latter approach within SPH 
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simulations dWhitehouse & Batel |2004 [2006; Whit ehouse et al.1 
120051: iMaver et al.ll2007[) and ap plied this both to large scale star 
formation simulations iBatd20091) an d to the fragmentati on of self- 
gravitating discs around young stars {Mem & Bate 201 1). There is 
however a substantial computational overhead associated with such 
calculations and this can impose undesirable limitations on the res- 
olution achievable 

An i ntermediate approach is r epresented by the al gorithm pro- 
posed by IStamatellos et al.l d2007h and developed by iForgan et al.l 
(2009). In its original form, this algorithm (which we henceforth re- 
fer to as 'polytropic cooling') computed a local cooling rate based 
purely on local variables]}} Specifically, it equates the cooling rate 
per unit mass with the local flux (F) divided by the column density 
(E) of material overlying a particular SPH particle. In the case of 
optically thin gas that is not subject to external irradiation, this re- 
duces to a cooling rate per unit mass that is simply ~ kctT a (where 
k is the opacity, cr the Stefan Boltzmann constant and T the tem- 
perature). This result is correct in the case of non-irradiated opti- 
cally thin gas in local thermodynamic equilibrium. In the case of 
optically thick gas in local thermdynamic equilibrium, the correct 
expression for the cooling rate per unit mass is V ■ F/p where p is 
the local density. As noted above, the polytropic cooling prescrip- 
tion instead sets this to F/E, which is not exactly the same but is 
likely to be of the same order provided that all quantities (including 
F) vary smoothly over a similar scale length. 

In order to calculate F/E from local variables, the polytropic 
cooling method makes a further assumption, that it is possible to es- 
timate E from the local value of the density and potential. In order 
to achieve this, it is assumed that the relationship between poten- 
tial Q¥), density (known for each particle in the SPH code) and E 
(to be estimated) can be represented by a mass weighted average 
over an equivalent n = 2 polytropic sphere in hydrostatic equilib- 
rium0 This estimator therefore purely uses the n = 2 polytrope as a 
smoothly varying spherical distribution from which one can conve- 
niently make a rough estimate of the relationship between p, *F and 
E. It does not require that the gas in the simulation is in hydrostatic 
equilibrium. In fact, the estimated relationship between p, *P and E 
(see Appendix) is very similar if one instead bases the mapping on 
a polytrope of different index, n, with the coefficient of E (Equation 
IA13t varying by aroun d 10% when n is varied over the range 0-5 
dStamatellos et alj20071) . Evidently, it is the smoothness of the dis- 
tribution, rather than the specific profile, that is the important prop- 
erty in this mapping procedure. Note that although the rationale for 
this approximation is based on the structure of smooth spherical 
distributions (i.e. where all quantities vary over a scale length of 
order the radius within the sphere) it has also been widely applied 
in simulations in which discs form and, as such, it is implicitly as- 
sumed that this mapping works in planar geometry. 

There are several variants of the above algorithm in use. For 
example, the cooling prescription also usually contains a term that 
prevents cooling below a background temperature (see Appendix). 
More recently, a hybrid method has been developed which com- 
bines the polytropic cooling method with flux limited diffusion 



1 Note that there is an obvious computational advantage in using quantities 
that are available for each particle and are pre-computed in the code. The 
addtional motivation for using the potential is that it contains information 
about the wider environment that may relate to the column density of over- 
lying gas 

2 The structure of this equivalent sphere is calculated only from the local 
density and potential and would therefore not in general be in hydrostatic 
equilibrium given the actual temperature of the particle. 



dForgan et alll2009h . This development follows from the recogni- 
tion that a possible drawback of the original method is that it al- 
ways implies cooling of gas (provided that it is hotter than the 
background temperature), whereas in optically thick regions, the 
sign of V ■ F might actually imply that a fluid element is heated by 
adjoining hotter regions. Flux limited diffusion treats this regime 
correctly but does not provide a good estimate of the cooling rate 
per unit mass in regions of low optical depth on the periphery of 
a condensation. It has been suggested that the advantages of both 
methods can be combined by simply adding the cooling rates from 
flux limited diffusion and from polytropic cooling. Clearly this can 
only be an improvement if each component of the cooling prescrip- 
tion makes a small contribution to the total cooling in the region in 
which this component is regarded as unreliable. This, however, has 
not been demonstrated to date. 

The development of such approximate cooling prescriptions 
represents an important opportunity to improve the verisimilitude 
of star formation simulations at a computational expense that is 
minimal compared with a full treatment of radiative transfer in the 
gas. It is obviously important that such prescriptions are rigorously 
tested. To date, the tests ha ve consisted of modelling th e collapse 
of a 1M spherical cloud dMasunaga & Inu tsuka 2000), th e sim- 
ulation of the collapse of a rotatin g cloud (|Bos s & Bodenheimer 
1979) and of a self-gravitating disc (Hubeny 1990]) and the thermal 
relaxation of a static spherical cloud (Spiegel 1957). In each case, 
the simulation results have been compared either with previous ra- 
diative transfer calculations or else with analytic results. Of these 
tests, only the latter, known as the Spiegel test, provides a direct 
measure of the cooling rates. This test (which involves the introduc- 
tion of sinusoidal temperature variations at a range of wavelengths) 
demonstrates that polytropic cooling provides a good measure of 
the cooling rate in the optically thin limit (as indeed it should, 
given that the cooling rate per unit mass in this limit is indepen- 
dent of the accuracy of the optical depth/column density estimate). 
However it cannot reproduce the results of this test (for arbitrary 
temperatute perturbation waevelengths) in the optically thick limit 
because in this case the cooling rate should depend on the wave- 
length, as this controls the rate at which radiation diffuses between 
hotter and colder regions. The polytropic cooling method, how- 
ever, gives cooling rates that are independent of the wavelength, 
since this rate is computed only using information about the local 
temperature and the estimated optical depth to the exterior, rather 
than on any temperature variations within the structure. The good 
m atch to the Spiegel test re sults in the optically thick limit found 
bv lStam atellos et al. ( 2007) in fact does not hold if one varies the 
wavelength from their adopted value. Note that this drawb ack has 
been addressed by the hybrid method o flForganetal]d2009l) (which 
correctly reproduces the results of the Spiegel test in the optically 
thick limit) since it is able to diffuse heat between hotter and cooler 
regions. 

The other two tests that have been conducted can instead be 
regarded as 'macroscopic' tests inasmuch as they look at the evo- 
lution of the entire system. A natural question that arises with such 
tests is whether their outcomes are as sensitively dependent on a 
correct treatment of the thermal physics as are the physical sit- 
uations that will be modelled with the algorithm. This is a hard 
question to answer unless one has also carried out 'microscopic' 
tests — in other words, an evaluation of the types of situation in 
which polytropic cooling does (or does not) provide a good esti- 
mator of the cooling rate. Once understanding of this question has 
been gained, one can then examine the physical situations that are 
encountered in real star formation simulations and decide whether 
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or not the algorithm is likely to be providing a reliable measure of 
the cooling rate. 

With this in mind, we here report on a suite of tests using sim- 
ple analytic structures in spherical and axisymmetric geometries. 
We first (Sj2} compare the predicted values of column density (£) 
and optical depth (t) provided by the polytropic cooling method 
with their true values obtained by integration of the analytic den- 
sity structures. We then proceed (ij3} to compare the polytropic es- 
timate of the cooling rate per unit mass (i.e. FJ1. where F s is the 
polytropic estimate of the radiative flux) with an equivalent quan- 
tity derived using true values of the column density and (in the op- 
tically thick limit) compare each of these quantities with what we 
term the 'true' cooling rate; this being related to the divergence of 
the radiative flux computed in the radiative diffusion approxima- 
tion. Section 4 summarises our conclusions. 

We caution at the outset that such a microscopic examination 
of the performance of the algorithm is designed to expose its weak- 
nesses and may expose weaknesses which actually turn out to be 
irrelevant in real simuations. If, for the sake of argument, we find 
a regime where the method produces cooling rates that are many 
orders of magnitude different from the true values, this may not ac- 
tually matter in practice. This is the case if, for example, both true 
and approximate values lie firmly in the regime where the gas is 
behaving adiabatically or in a regime where the gas temperature 
is controlled by the ambient cloud temperature. Here we simply 
present our results in order that those who use this method in hy- 
drodynamical simulations can assess the significance in practice of 
the discrepancies that we highlight here. 



2 COMPARISON OF COLUMN DENSITY AND 
OPTICAL DEPTH ESTIMATES 

2.1 The application of the polytropic method to spherical 
hydrostatic structures 

We set up a set of simple model density and temperature distri- 
butions corresponding to various polytropes in hydrostatic equilib- 
rium. Hydrostatic structures are of interest when testing the radia- 
tive cooling approximation since these systems will, once formed, 
tend to exist for prolonged periods in hydrodynamical simulations 
and their thermal evolution in a given state will be important as well 
as their dynamical evolution. 

Given an assumed dependence of the Rosseland mean opac- 
ity on density and temperature, we can readily compute the col- 
umn density (£) and optical depth (t) along a radial path from any 
given point to infinity. We then compare these quantities t the col- 
umn densities and optical depths predicted by the polytropic cool- 
ing method (L,, i y and T po i y ; see Appendix). We emphasise that T pD / Y 
is not simply the product of the local opacity and the estimated 
column density but also contains a (opacity law dependent) factor 
that estimates the ratio of local opacity to the pat h averaged opacity 
within the polytropic smoothing formulation; see lStamatellos et all 
(2007) and Appendix IA4I This aspect of the polytropic cooling for- 
mulation is designed to improve the optical depth estimate in cases 
where the opacity is a strongly varying function of temperature 
(and hence position) as in the so-called 'opacity gap' in protostellar 
discs. 

It is worth noting some features of the approximation that will 
be useful in interpreting our results. Members of a polytropic fam- 
ily of given n can all be mapped onto a single function representing 
the dependence of a scaled density-like variable on a scaled radius 




Figure 1. Density profile and cumulative mass distribution for the n = 2 
polytrope (in arbitrary units) which is the basis for the mapping between p, 
¥ and X in the polytropic cooling method. Note that the mass distribution 
is weighted towards larger radius. 



variable (<f). The polytropic cooling approximation first assumes 
that a given particle corresponds to a particular value of £ within 
a polytropic sphere and then, given the actual values of density 
and potential at that particle, determines the corresponding value 
that the column density would have in that case. This is repeated 
over a range of <f values and then a final 'effective' column den- 
sity is computed which weights the contributions from different £ 
in proportion to the amount of mass at different £ values within the 
polytropic sphere. Since for an n = 2 polytrope, the majority of 
the mass is located at rather large radii (see FigureQ}, this implies 
that the relationship between f, p and £ is equivalent to assum- 
ing that the particle resides towards the outskirts of a structure. It 
therefore implies a lower value of £ (for fixed *F and p) than for 
a particle that was mapped onto the centre of a polytropic sphere. 
The polytropic cooling method applies this weighting over f to all 
particles since, being entirely local, it has no knowledge of where 
any particular particle is located within a parent structure. There- 
fore if, in fact, a particle is close to the core of its parent structure, 
the weighting of the polytropic cooling method will systematically 
underestimate the overlying column density. Conversely, for parti- 
cles in the extreme periphery of a structure, the weighting of the 
polytropic method over-estimates the column density at given f 
and p. 

This simple insight is sufficient for us to understand the 
method's systematic under-estimate of the column density near 
the centre of spherical structures and, conversely, its systematic 
over-estimate in the outer regions (Figure The method under- 
estimates the central column density by about a factor three for all 
polytropic spherical structures and the over-estimate at large radius 
is around a factor three for polytropes with n of 2 or greater and 
larger for smaller n (less compressible gas). 

In the case of constant opacity (e.g. that due to electron scat- 
tering), these results for the column density apply directly to the 
derived optical depth since the two quantities are simply propor- 
tional. In Figure [3] we present the optical depth as a function of 
radius for the case of opacity given by ice grains for which k oc T 2 . 
Note that normalisation of optical depth is arbitrary. This illustrates 
that within the region containing the bulk of the mass, the poly- 
tropic method underestimates the optical depth. 
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Figure 2. Algorithm tests for spherical polytropes. Upper row: n = (uniform density), middle row: n = 1 polytrope, lower row: n = 5 (nearly isothermal) 
polytrope. In each case, the left hand panel is the density profile, the middle panel compares the integrated column density to infinity from each point (solid 
line) with the corresponding quantity derived by the polytropic cooling method (dashed line) and the right hand panel is the ratio of values for the column 
density on a logarithmic scale. Units are arbitrary. 




Analytic 

Polytropic method 
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Figure 3. Comparison of the optical depth (in arbitrary units) due to ice 
and the polytropic estimate of the same quantity (dashed line) in the case 
of the n = 4 polytrope. Note that the normalisation for the optical depth is 
arbitrary. 



2.2 The application of the polytropic method to disc-like 
geometries 

We test the polytropic method in planar geometry, being mindful 
of the fact that this applies to a number of situations that occur 
in hydrodynamical simulations such as shock compressed sheets 
and rotationally supported discs0 To this end we use an approx- 
imately Keplerian disc model in which both the column density 
(X) and vertically averaged temperature vary as R 3 . Such a struc- 
ture h as a radially constant profile of Toomre Q parameter jToomrej 
1 1964b as would be expected in the case of discs that are in a state of 
self -regulated margina l stability against self-gravity dGammie et al.l 
120001 : iRice et aLl botBl. In fact, this particular profile corresponds 
to the case of a disc with constant accretion rate and rad iative cool- 
ing with opacity provided by ice grains I Clarke! 2009), an opac- 
ity regime that is applicable to the outer regions of discs around 

3 Note that although Stamatellos & Whitworth (2008) demonstrated that 
the method satisfies the Hubeny test ? for the variation of temperature as a 
function of optical depth within a disc, this does not demonstrate that the 
method has calculated the optical depth correctly, which is the focus of our 
investigation here. 
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young stars and where the cooling rate is important in determining 
whether such discs fragment. 

Although (as emphasised by Stamatellos et al 2007), their 
cooling prescription is not only applicable to situations that are 
in hydrostatic equilibrium, we here test the method on a density 
and temperature field that most closely resembles the situations that 
arise in star formation simulations in which hydrostatic equilibrium 
normal to the disc plane is attained on a short (dynamical) timescale 
(we assume that the small influence of pressure gradients and of the 
disc's self-gravity in the disc plane is balanced by a minor adjust- 
ment of the disc's rotation curve). In order to obtain an analytic 
approximation to the disc vertical structure as a function of radius, 
we approximate the vertical gravitational acceleration at height z in 
the disc as 



gz 



-n 2 z - AtxGI.' 



(i) 



where the first term is the vertical component of the gravitational 
acceleration provided by the central star (CI is the local Keplerian 
angular velocity) and S' is the column density between the point 
at height z and the disc mid-plane (which is not to be confused 
with the column density to infinity that has been discussed hith- 
erto). This second term is the gravitational acceleration that the 
disc would exert in the limit that it was an infinitely extended strat- 
ified slab. In reality there are further small corrections to the self- 
gr avity due to radial g radients in the disc (see e.g. the appendix 
of iBertin & Lodatolll999h . We have verified a posteriori (by sub- 
tracting the density distribution of a radially constant disc from the 
computed structure and then calculating g z by direct summation) 
that this further correction would account for at most ~ 3% of the 
total g z value. We emphasise that this approximate g z value is only 
used in order to obtain a structure that is in approximate vertical 
hydrostatic equilibrium and that the full disc (and star) potential is 
used when computing the polytropic cooling. 

The motivation for adopting this approximate form of g z is 
that it allows us to compute an analytic disc structure as a function 
of cylindrical radius (R) and z, since such a disc can be modeled by 
an n = 2 polytrope with P = Kp 2 where K is a global constant^ In 
this case, differentiation with respect to z of the hydrostatic equilib- 
rium equation normal to the mid-plane yields the modified simple 
harmonic equation: 



dz 2 



2nG, Q. 2 



whose solution (subject to p = p and p' = at z = 0) is 
P=Pof(l+/(Q))cosJ--/(Q) 
where 



/(G) : 



and 



n 2 



AitGpa 



H0 = {^) 



1/2 



(2) 



(3) 



(4) 



(5) 



4 Note that although this equation of state yields the surface density profile 
of a steady state marginally self-gravitating disc with opacity dominated 
by ice cooling, it does not do so uniquely, since the actual pressure-density 
relationship in such a disc depends on issues such as the vertical profile 
of viscous dissipation in the disc. It nevertheless offers a simple analytic 
desription of a disc with an astrophysically relevant surface density profile. 
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Figure 4. The column density in the z direction out to infinity in disc geom- 
etry, compared to the polytropic estimate. 



The constant / is of order Q when Q » 1 (the non-self grav- 
itating case) while / ~ Q 2 in the strongly self-gravitating case 
(Q ^ !)■ In the marginally self-gravitating (self-regulated) case 
that we consider here, / is of order unity. 

The density profile given above allows us to derive the scale 
height of the disc (defined here as being the height at which the 
density falls to zero): 



H = Hqcos 



AQ) ) 
1+/(G)J 



(6) 



We set up such a disc over the radial range R = 1 to R = 10 (in 
arbitrary units), and adopt parameters / = 1 and H = 0.35. For this 
model the total disc mass (which, given the surface density profile, 
is concentrated at small radii) is ~ 20% of the central star mass. 
At given R and z, we compute the gravitational potential by adding 
the contributions due to the star and the rest of the disc. The latter 
is calculated by using elliptic integrals to compute the contribution 
to the potential from each ring of material corresponding to each 
pair of R and z values in the grid. The local density and potential is 
then used to obtain the polytropic cooling method's estimator of the 
column density which is compared with the 'true' column density 
integrated out to infinity in the z direction. Figure[4]compares these 
column density values as a function of z for a cylindrical slice at 
R = 5, at which point the disc aspect ratio H/R ~ 0.07. The upper 
curve in the Figure represents the case that the full (star plus disc) 
potential is used in the estimation of E, whereas the intermediate 
curve uses just the disc potential (an approach suggested by Sta- 
matellos et al 2007 though not always implemented in disc studies). 
As expected, the estimated column density is larger when the stel- 
lar potential is included by around a factor 2 ~ 3, consistent with 
the square root dependence of the estimated column density on the 
potential and the fact that the disc mass is concentrated towards the 
inner part of the disc. However, even in the case that only the disc 
potential is included, the estimated column density exceeds the true 
column density by a factor 5 or more at all heights. The consistent 
sign of the discrepancy (in the sense that the estimated column is 
always larger) is simply because in a disc geometry there is always 
a low column density route to infinity (along the symmetry axis) 
which is not accounted for by the estimated value. 
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2.3 Summary 



3.1 Cooling rates in spherical geometry 



We find that in the case of spherical structures, the polytropic 
method under-estimates the column density (and hence optical 
depth in the case of constant opacity discs) by a factor of 3 ~ 4 
in the inner regions and over-estimates the column density at large 
radius. This can be understood in terms of the weighting function 
that is used in order to map the density and potential onto the col- 
umn density which is weighted towards the relationship between 
these quantities that applies at intermediate radii within an n = 2 
poly trope. 

On the other hand, the method systematically over-estimates 
the (vertical) column density in a disc structure, even in the case 
that one uses only the potential from the disc (rather than the disc 
plus star) in estimating the column density. For the disc structure we 
have modeled (for which H/R ~ 0.07 at the radius considered and 
where Q T ~ 1 throughout), this over-estimate is around a factor 
5 or more. This over-estimate can be understood as a geometric 
effect in that the estimate based on local density and potential does 
not take account of the direction of low £ normal to the disc plane. 



3 IMPLICATIONS OF RESULTS FOR THE ESTIMATE 
OF COOLING RATES 

The polytropic method uses the column density and optical depth 
estimates discussed above to compute the cooling rate per unit mass 
according to 



U = 



crT 4 



(7) 



Here, r is the path averaged optical depth (see Appendix IA3t 
whereas f is the product of local opacity and estimated column 
density. These two quantities differ by a factor that is usually of 
order unity, except in regions where the opacity is strongly depen- 
dent on temperature (see Figure 6 of lStamatellos et al. 2007). Note 
that in both cases the opacity is given by parametrisations of the 
Rosseland mean opacity 

There are several approximations involved in Equation[7j First 
of all, it is assumed that the cooling rate per unit mass at a given 
point is simply an estimator of the local radiative flux divided by the 
estimated column of material above this point. Secondly, it takes a 
form for the local radiative flux designed to have particular func- 
tional forms in the limits of small and large optical depths. At low 
optical depths this is simply the product of the optical depth and the 
black body flux. In the limit of large optical depth, this expression 
is motivated by the radiative diffusion approximation 

4 



3Kp 



v(o-r') 



(8) 



and is likely to be of similar order to Equation [8] provided that all 
quantities vary over a scale length r (so that — ~ ^p). Finally, 
the second term in the numerator of Equation [JJ ensures that the 
cooling rate goes to zero at temperature T and is designed to mimic 
the effect of an external radiation field at this temperature. For the 
purposes of the test we set T = 0. 

However, it is important to recognise that Equation [7J is itself 
an approximation. In the optically thick limit we can also evaluate 
the 'true' cooling rate per unit mass, using 



1 

-V-F 

P 



(9) 



where F is evaluated in the diffusion limit (i.e. Equation[8). 



We first examine how the cooling values predicted by Equation 
[7j using column density and optical depths estimated by the poly- 
tropic method compares with the same quantity evaluated using the 
actual optical depths and column densities. In each figure, the left 
hand panel refers to the case of constant opacity and the right hand 
panel to the case k oc T 2 . For the purpose of this comparison, we 
evaluate the temperature profile (Figure [5} by assuming that the 
structure is in hydrostatic equilibrium. The (true) optical depth to 
the centre is 42 for the constant opacity case and 100 for the case of 
ice and the point at which the (true) optical depth is unity is marked 
by the arrow in Figure [3~T1 

Figure [3~T| shows that the polytropic cooling estimate of Equa- 
tion[7J(black line) has the same general shape as the same quantity 
evaluated from the actual column density (green dotted line), in that 
in both cases the cooling rate per unit mass is highest close to an 
optical depth of unity and that both converge to the same value in 
the limit of low optical depth. This is to be expected since as t tends 
to zero, the cooling rate per unit mass tends to a value that is inde- 
pendent of column density and depends only on temperature (the 
optically thin limit of Equation [7J is shown by the blue dash-dot 
line in Figure [3~lt . However, in the optically thick portion of the 
disc, the two cooling rates are in general quite different, with the 
polytropic estimate over-estimating the cooling at small radii and 
under-estimating it at large radii. This can be simply understood in 
that the cooling rate per unit mass in the optically thick limit scales 
as 1 /E 2 and hence this result reflects the respective under-estimate 
and over-estimate of £ at small and large radius. 

So far, we have just examined how the accuracy of estimating 
the column density affects the cooling rate as specified by Equa- 
tion[7J A more important comparison is however the comparison of 
the polytropic method's estimated cooling rate with the 'true' (ra- 
diative diffusion) method in the limit of high optical depth (Equa- 
tions [8] and [9}, shown as the dashed red line in Figure [3~T| Inter- 
estingly, this agrees rather well with the polytropic method in the 
optically thick region (the black line) in both opacity regimes. In 
both cases the total (mass integrated) cooling rates agree to within 
a few tens of per cent. Such a discrepancy is almost certainly no 
larger than additional errors due to poorly determined opacities, for 
example. We have also explored a variety of different polytropic 
spherical structures and different opacity laws and find that the 
polytropic cooling method's fair (i.e. order unity) agreement with 
the results of radiative diffusion is reproduced for a wide range of 
input parameters. This explains why the method performed well 
in reproducing the results o f the spherical collapse calculation of 
iMasunaga & Inutsukal ( feOOOh . 

When the cooling rates are integrated over mass through the 
optically thick part of the structure, one obtains the result that the 
polytropic method modestly under-estimates the total cooling rate 
by around 20% (compared to the radiative diffusion limit) in both 
cases. The switch in the sign of the discrepancy at intermediate 
radius means that the global values agree rather well despite the 
discrepancies in local rates. 

It is perhaps surprising that the polytropic cooling method re- 
produces the 'true' (radiative diffusion) cooling rate considerably 
better than it reproduces the quantity that it ostensibly computes 
(Equation IT) . This is because Equation |7j is only a reasonable ap- 
proximation of Equation [8] in the case that one can replace differ- 
entiation with respect to t by simple division by r. Comparison 
between the green and red lines in Figure [37X1 show that this is a 
poor approximation at small radius, i.e. the temperature falls over 
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Figure 5. Density (left panel) and temperature (right panel) profiles, in arbitrary units, of the hydrostatic n = 4 polytrope used for the cooling rate comparisons. 
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Figure 6. Comparison of estimates of the cooling rate per unit mass (in arbitrary units) for the case of a) constant opacity and b) ice grain opacity for 
n = 4 polytropic spheres with (true) central optical depths of 42 and 100, respectively. The black line represents the polytropic estimate to Equation 
UJ while the green dotted line uses the same equation but with the true values of the optical depth and column density. The blue dash-dot line is the 
limiting form of both expressions at low optical depth where they agree by construction. The red dashed line is the 'true' cooling rate computed using 
the radiative diffusion approximation (Equations [S]and|9j. The arrows denote the point of (true) optical depth equal to unity. 



an optical depth that is a small fraction of the total optical depth 
at small radius and hence simple division by t under-estimates the 
magnitude of the cooling. This is almost exactly compensated by 
the fact that the polytropic cooling method under-estimates the op- 
tical depth in the central regions (see ij2j. Although it may be re- 
garded as unfortunate that the reason for this good agreement is 
a near cancellation of discrepancies of comparable magnitude, we 
are encouraged that, as described above, the reasonable agreement 
appears to hold - in spherical geometry for a wide range of input 
parameters. 

3.2 Cooling rates in disc geometry 

The situation in disc geometries is very different. As above, we 
start by comparing the polytropic cooling estimate of Equation [7J 
with the same quantity evaluated using the true column density at 
each point. For consistency with the cooling regime encountered in 
the outer regions of self-gravitating discs, we employ the opacity 
appropriate to ice cooling (k oc T 2 ). We here focus on the optically 
thick regions of the disc, in order to allow comparison with Equa- 
tion [8] and so neglect the additive term oc f" 1 in the denominator 
of Equation[7] Figure [TJplots the specific cooling rate estimates as 



a function of z for the slice at if = 5 in the model disc described 
in 92.21 The dashed line represents the cooling rate according to 
the radiative diffusion approximation (i.e. Equations [8] and |9) while 
the upper solid line employs Equation [7J where the column den- 
sity is the 'real' value and the optical depth is approximated as 
the product of the real column density and the local opacity. The 
lower two lines correspond to Equation fTJ evaluated according to 
the polytropic method where the optical depth is the product of the 
estimated column density and the estimated average opacity along 
the line of sight (see Appendix I A4t . In the case k oc T 2 this is 0.58 
times the local opacity. The bottom line uses the full (star plus disc) 
potential in estimating the column density whereas the curve above 
uses the disc potential only. 

Recalling that in the optically thick regime the, cooling rate 
according to Equation [7] scales as E~ 2 and that the method over- 
estimates the column density in disc geometries, it is no surprise 
that now the polytropic cooling method yields a much smaller cool- 
ing rate, being at least an order of magnitude lower in the case that 
only the disc potential is applied and more than two orders of mag- 
nitude if the total (star plus disc) potential is used. 

When we compare each of these rates with the 'true' cooling 
rate in the radiative diffusion approximation (dashed curve, Equa- 
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Figure 7. The cooling rate for the disc system evaluated by the polytropie 
cooling method (both including and excluding the stellar potential) com- 
pared with the diffusion limit as well as Equation [JJ using analytic optical 
depths at radius R = 15. 

tionO we find that the latter is now of comparable magnitude to the 
cooling computed from Equation [7] using the true column density 
but is far larger than the estimates from polytropie cooling. Thus in 
contrast to the case in spherical geometry, there is no compensating 
effect that brings the polytropie estimate into agreement with the 
'true' cooling rate. 

Figure[8]illustrates how the specific cooling rates that are com- 
pared in Figure [7] affect the total cooling rate per unit area. At each 
value of z, the quantity plotted is given by the integral of the prod- 
uct of specific cooling rate and density from z' = to z. Note that 
in contrast to Figure [7] Figure [8] is plotted on a linear scale. It il- 
lustrates that at all heights, the use of Equation [7] gives a reason- 
able (order unity) estimate of the integrated cooling rate provided 
one has knowledge of the 'true' column density. Indeed, the total 
cooling rate per unit area is within 25% of the result of using the 
radiative diffusion approximation. The total cooling rate from the 
polytropie method is, however, too low by a factor of at least ~ 200. 
This emphasises that, in disc geometry, the most pressing problem 
is to find a reliable estimator of the column density. 

We have also experimented with varying the opacity law in the 
disc geometry and note that when the opacity is a less strongly in- 
creasing function of temperature, there is a greater tendency (com- 
pared with the spherical case) for discs to be in the regime of net 
heating (i.e. the upper regions receive a greater flux from the warm 
disc below than they radiate upwards — the outwardly increasing 
area element in spherical geometry makes this less likely in the 
spherical case). Obviously, no local cooling prescription can repro- 
duce this effect. Thus any future efforts to find an improved local 
cooling prescription in disc geometries will need to be carefully 
calibrated in different opacity regimes. 



4 CONCLUSIONS 

The polytropie method works quite well in spherical geometry. We 
have shown that over a wide range of input parameters the mass in- 
tegrated cooling rate out to the r = 1 surface matches that from ra- 
diative diffusion to within order unity or better, an error that is prob- 
ably negligible compared with other uncertainties such as those in 
the opacities. This good agreement is 'fortuitous' in the sense that 
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Figure 8. Magnitude of the integrated cooling rate per unit area for the disc 
evaluated by the polytropie cooling method compared with the diffusion 
limit as well as Equation[7Jusing analytic optical depths at radius R = 15. 



it results from two nearly compensating discrepancies, but can nev- 
ertheless be exploited in simulations in spherical symmetry. 

The polytropie cooling estimate, however, performs poorly in 
optically thick discs, regardless of whether one uses the potential 
from the disc alone or that due to the star and the disc in order to es- 
timate the column density. This is due to the fact that the polytropie 
method over-estimates the column density considerably (ijjjl since 
it does not take account of the lower column density route normal 
to the disc plane. This translates into an under-estimate of the cool- 
ing rate by between two and three orders of magnitude (Figures [7] 
and[H), depending on whether the stellar potential is included in 
the column density estimate. We have shown however that (at least 
in the astrophysically relevant ice opacity regime) the cooling can 
be well modelled in disc geometries if only one can improve the 
column density estimate. This should clearly motivate further algo- 
rithmic developments in this area. 

The above is sufficient to say that the method needs to be used 
with caution, though it does not necessarily negate results in the 
literature that use this method. For example, the outer regions of 
proto-stellar discs may in fact be nearly isothermal at the back- 
ground temperature, in which case all methods would agree on the 
thermal structure of the disc. Alternatively, even if simulations en- 
ter regimes where the cooling rate is strongly under-estimated, it 
may have a negligible effect on where, for example, fragmentation 
occurs in the disc. This is becaus e the cooling rate in such discs is 
a very strong function of radius ( Clarke 2009) and thus an under- 
estimate in cooling would cause only a modest increase in the ra- 
dius at which the disc fragments. 

Simulators evidently need to balance these issues against the 
computational economies offered by polytropie cooling. These re- 
sults can also inform those wishing to combine polytropie cooling 
with other methods, for example the expedient of adding the poly- 
tropie cooling method to flux limited diffusion (as in the 'hybrid' 
method of Forgan et al) relies on the former being a minority con- 
tributor in the optically thick regions where the latter is correct. 
We find that this requirement is amply fulfilled. On the other hand, 
whether or not the hybrid is an improvement over pure flux limited 
diffusion depends on the situation. Unlike the flux limited diffusion 
method, the hybrid method would yield the correct cooling rate for 
non-irradiated optically thin gas. Neither method, however, is suit- 
able for treating the optically thin regions above an optically thick 
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disc which are subject to irradiation from the disc below. Whether 
or not this matters for the system dynamics obviously depends on 
the topic at hand. 

Finally, we draw attention to the fact that an astrophysically 
relevant disc structure (a marginally self-gravitating optically thick 
steady state disc with opacity dominated by ice grains) can be mod- 
elled as an n = 2 polytrope whose density and temperature structure 
can be described analytically. The existence of such an analytic so- 
lution may be useful in testing future cooling algorithms. 
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APPENDIX A: THE POLYTROPIC COOLING METHOD 

The radiative cooling of an element within a cloud is determined by 
the optical depth of that point, that is the integrated opacity along 



th e line-of-sight out of the cloud. The polytropic cooling method 
of IStamatellos et al] d2007l) aims to evaluate the radiative cooling 
of an element by calculating an approximate optical depth using 
only local variables at the point of interest. 

In order to achieve this, the particle of interest is taken to 
be embedded at an arbitrary position in a spherically symmetric 
pseudo-cloud, modelled by a polytrope of index n. The density (p), 
potential (f) and temperature (T) of the polytrope are given by the 
solutions of the Lane Emden equation, 9„, 



-47TGR 2 oPc m 



(Al) 
(A2) 
(A3) 



where = -&(§) p=& + 0(f). 

The parameters of the polytrope are set to reproduce the local 
conditions at the location of the particle. The column density, opti- 
cal depth and radiative cooling are then calculated for the particle 
embedded in the pseudo-cloud. The value of the column density to 
the particle of interest is taken to be the average column density 
over all possible positions within the pseudo-cloud. 



Al Calibrating the Pseudo-Cloud 

Once the polytropic index of the barometric equation of state is set, 
the values of the central density, scale length and core temperature 
for the pseudo-cloud in which the particle of interested is embedded 
at a radius r = Rq^ are calculated such that the known density (p,), 
gravitational potential pP,-) and temperature (T,) at the position of 
the particle are reproduced. 



Pc 

Ro 



(A4) 
(A5) 
(A6) 



A2 Column Density 

The column density (mass per unit area along a line-of-sight) from 
a particle at radius r = £Ro along a radial line to the edge of the 
pseudo-cloud = (; B ) is given by 



2 = J " p(.r)dr = R J p(f)df 
Substituting for the density using Equation lAll 



2(f) = i? J PcPigW 



And substituting for p c and R using Equations lA4l and |A5| 



(Al) 



(A8) 



(A9) 



When calculating the mean column density over all positions 
of the particle within the pseudo-cloud, the column density at a 
given radius E(f ) is weighted by the mass contained in a spherical 
shell radius £, thickness df , M(f ) as this gives the probability of the 
particle being within that shell. The total mass of the cloud, M is 
used to normalise the distribution. 



1 f f8 
M Jo 



m)M(t)cU; 



(A10) 
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in 


0.0 


0.3796 


1.0 


0.3756 


2.0 


0.3682 


3.0 


0.3600 


4.0 


0.3500 


5.0 


0.3322 



Table Al. Values of f„ for different polytropic indices (evaluated numeri- 
cally). 



The mass enclosed in a spherical shell is given by the volume of 
the shell multiplied by the density, from the Lane-Emden func- 
tion, M(g) = 47tRlp c fe"(^). The total mass of the cloud is M = 



Substituting these in to the integral gives 



Jo ^ 



(g)dg 



(All) 



Evaluating the first integral and substituting in the column density, 

-l 



£ = 



-e B 



Which can be written 



Xf« r r(B 
J ng)dg 



-W(f 

AnGm 



£ = 



-%Pi 
4nG 



f 2 ^(A12) 



(A 13) 



Where is just a constant for a given polytropic index 

-i-i 



, 2 d0 



Mb 



rs r r^B 



0(f) 



Values of can be evaluated numerically and are given in Table 
IA2I It can be seen that the results for the column density are fairly 
insensitive to the polytropic index chosen in the approximation. 

This mean is used as the approximate value of the column den- 
sity at the point in question and depends only on the local potential 
and density. In this formula, the local potential conveys the struc- 
ture of the surrounding space (rather than performing the integral 
along the line-of-sight). 



The mass-weighted pseudo-mean optical depth (averaging 
over all possible positions of the particle within the pseudo-cloud) 
can be obtained as for the column density. 



-f. 






If) 






r=f«J 


-(. 






I?) 





-%pi 



4nG 



J^b r /^s / 



0(g) 



e(g) 



e"(g)dg 



A4 Pseudo-Mean Mass Opacity 

Once the column density and optical depth have been calculated, 
one can define an average opacity along the line-of-sight, 

f 

K = — 



k(p,T) 



-Ub 



«R P 



0(f) 



0(g) 



0(g) 



0(0 



new 



Q"(g) 



k is a function only of p and T, the local density and temperature, 
and does not depend on particle positions, therefore it need not be 
evaluated on-the-fly during simulations but can be evaluated in ad- 
vance and looked-up from a table during simulations. Multiplying 
by the pseudo-mean column density, E, recovers the optical depth. 



A3 Optical Depth 



The optical depth along a ray path with absorption coefficient a„ is 
defined as t„ = j a u (r)dr with a u = K u p. This can be averaged over 
all frequencies in the black body spectrum using the Rosseland- 
mean opacity Kg(p, T). 

The optical depth of a particle at radius r = ^R in the pseudo- 
cloud can be calculated, 



r 



r(g) = K R (p(g),T(g))p(g)R dg 



(A14) 



Substituting for the density and temperature using Equations IA1I 
andlA3l 

T(f) = J K R ( Pc 0'(g\ r c 0(f )) Pc 0"(f )# <2f 



(A15) 

And substituting for p c , T c and R using Equations lA4l|A5| and |A6| 

9(g)" 



T(f> = 



4nG<KO 



0"(g) ( 



0(g) 



-T, 



0(g) 



0(g) 



0(0 



dg (A16) 
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